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ABSTRACT 


A  model  for  ELF  atmospheric  noise  is  proposed  and  studied  in  detail.  The 
model  is  completely  described  by  a  few  meaningful  parameters  which  are 
easily  correlated  with  measured  data.  Procedures  for  fitting  the  free  param¬ 
eters  to  observed  noise  data  are  provided.  The  model  should  prove  useful  in 
the  analysis  and  design  of  communication  systems  operating  in  impulsive 
noise  environments. 
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A  MODEL  FOR  ELF  NOISE 


I.  INTRODUCTION 

The  primary  source  of  electromagnetic  noise  at  ELF  (3  to  300  Hz)  tan  be  attributed  to 
lightning  discharges  which  propagate  in  the  earth- ionosphere  cavity.  As  a  result,  the  received 
noise  at  a  typical  ELF  receiving  site  tends  to  be  impulsive  in  character  even  in  the  absence  of 
strong  local  thunderstorm  activity.  Furthermore,  ELF  noise  has  been  observed  to  exhibit 
strong  geographical,  seasonal,  and  diurnal  dependence  so  that  the  observed  noise  background 
can  be  expected  to  be  non- stationary  as  well  as  highly  non-Gaussian.  To  assess  the  perform¬ 
ance  of  communication  systems  operating  in  such  environments,  it  is  important  to  have  available 
a  mathematical  model  that  provides  a  complete  statistical  description  of  the  received  noise.  It 
might  be  added  that  this  requirement  exists  whether  the  problem  is  that  of  obtaining  results  on 
system  performance  analytically  or  by  computer  simulation.  To  be  useful,  any  such  model  must 
satisfy  the  following  discrata: 

(a)  It  must  agree  reasonably  well  with  the  underlying  physics. 

(b)  It  must  be  completely  describable  by  a  few  meaningful  parameters. 

(c)  Its  statistical  characterization  must  be  easily  correlated  with  corre¬ 
sponding  measured  data. 

(d)  Hopefully  it  i3  mathematically  tractable. 

The  present  report  is  one  of  a  series  concerned  with  the  development  and  application  of  such 
a  model  to  ELF  communications.  In  this  report  we  will  be  concerned  with  the  development  of  this 
model  and  the  fitting  of  its  parameters  to  observed  ELF  noise.  Other  reports  in  this  series  will 
be  concerned  with  optimum/suboptimum  receiver  structures  and  their  digital  realizations. 

Before  entering  upon  a  detailed  discussion  of  this  model,  let  us  briefly  review  some  of  the 
pertinent  characteristics  of  ELF  noise. 

H.  CHARACTERISTICS  OF  ELF  NOISE 

General  observations  of  wideband  ELF  noise  reveal  frequent  occurrences  of  large  pulses 
superimposed  upon  a  more  homogeneous  appearing  background  noise.  Typical  recorded*  H-field 
noise  waveforms  are  illustrated  in  Figs.  1  through  3.  We  have  found  it  convenient  to  classify  the 
recorded  noise  into  the  three  categories:  high,  moderate,  and  low-level,  depending  upon  the 
observed  spectral  levels.  An  exception  is  the  case  of  Norway  data  which  appear  to  be  of  a  uni¬ 
form  level.  The  large  pulses  present  in  the  waveforms  of  Figs,  i  through  3  are  undoubtedly  due 
to  strong  local  thunderstorm  activity  while  the  background  structure,  on  the  other  hand,  is  due 
fo  a  considerably  greater  number  of  unresolved  low- level  pulses  attributed  to  more  remote 
thunderstorm  activity.  As  a  result,  it  follows  from  central  limit  theorem  considerations  that 
the  large  number  of  low-level  overlapping  pulses  resemble  Gaussian  noise.  Considerable  de¬ 
partures  from  this  condition  are  to  be  expected  at  larger  amplitudes.  Indeed,  measurements 
of  the  amplitude  probability  distribution  (APD)  of  ELF  noise  confirm  this.  Figures  4  through  6 
illustrate  the  computed  APD's  corresponding  to  the  noise  waveforms  appearing  in  Figs.  1 
through  3.  Here  we  are  plotting  the  probability  or  percent  of  time  that  the  absolute  value  of  the 
recorded  noise  exceeds  a  given  level  measured  in  rms  units.  The  scales  have  been  chosen  so 
that  the  APD  of  Gaussian  noise  plots  as  a  straight  line  of  slope  as  illustrated  in  Fig.  V,  From 

*The  details  of  the  measurement  and  recording  system  are  described  In  Ref.  1. 
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Fig.  3.  Typical  H-field  ELF  noise  waveform  recorded  at  Norway. 

Figs.  4  through  6  it  is  to  be  observed  that  at  low  amplitudes  the  APD  does  indeed  resemble  that 
of  Gaussian  noise  while  exhibiting  a  much  larger  dynamic  range  at  larger  amplitude  levels. 

Also  of  interest  are  various  first-order  temporal  statistics  of  ELF  noise.  We  will  be  par¬ 
ticularly  interested  in  such  quantities  as  the  average  number  of  high-level  pulses  per  second  and 
the  distribution  of  the  interval  between  their  occurrences.  Unfortunately,  such  quantities  are 
difficult  to  obtain  from  the  recorded  data  in  any  unambiguous  fashion.  Some  techniques  for  the 
analysis  and  measurement  of  the  temporal  statistics  of  ELF  noise  are  described  in  detail  in 
Ref.  i.  Here  a  useful  criterion  for  the  presence  of  a  pulse  in  wideband  ELF  noise  was  presented 
in  terms  of  three  parameters  TL,  Tu,  and  Sc>  as  illustrated  in  Fig.  8.  The  "beginning"  of  a 
pulse  is  said  to  have  occurred  when  the  absolute  value  of  the  noise  waveform  first  exceeds  a 
threshold  T  ,  while  the  "duration"  of  such  a  pulse  is  measured  from  this  instant  until  such  time 
as  the  absolute  value  has  been  below  for  at  least  Sc  seconds.  It  is  clear  that  any  practical 
usp  of  measurements  based  upon  this  criterion  will  require  careful  choice  of  the  parameters 
T^,  Tu,  and  Sc»  Obviously,  both  and  Tu  should  be  placed  just  above  the  background  noise 
level  if  such  a  point  can  be  accurately  defined.  Often  this  can  be  accomplished  by  observation 
of  the  APD  plot.  The  knee  of  the  APD  plot,  which  represents  the  point  of  departure  from 
Gaussian  behavior,  is  in  most  cases  quite  pronounced  (Figs.  4  through  6).  Suitable  thresholds 
taken  about  this  point  should  provide  useful  information  on  the  temporal  statistics  of  the  high- 
level  pulses  while  excluding  the  low-level  background  noise  from  consideration.  For  definiteness 
we  have  developed  the  following  useful  procedure  for  placement  of  the  threshold  levels.  First, 
draw  tangents  to  the  two  distinct  portions  of  the  APD  curve  as  illustrated  by  the  dashed  lines  in 
Fig.  9.  Let  T  represent  the  level  in  rms  units,  corresponding  to  the  intersection  of  these  two 
lines  (T  =  -  6  dB  for  the  high-level  Saipan  data  illustrated  in  Fig.  9)  and  set  T  =  T.  =  T.  The 
choice  of  Sc,  on  the  other  hand,  is  not  quite  so  straightforward.  Obviously,  we  would  like  to 
make  Sc  as  small  as  possible  without  causing  portions  of  a  single  noise  burst  to  be  mistaken 
for  two  or  more  resolved  pulses.  From  experimentation,  with  the  noise  thresholds  set  at  the 
value  T  described  above,  we  have  found  values  of  Sc  in  the  range  5  to  10  msec  sufficient 
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TABLE  I 

TYPICAL  FIRST-ORDER  TEMPORAL  STATISTICS  OF  ELF  NOISE 


Average  Number  of  Pulses/Second 

Average  Pulse  Duration* 

Noise  Data 

Sc  =  5  msec 

S  =  1 0  msec 
c 

=  5  msec 
c 

S_  =  10  msec 
c 

Saipan 

High-  Level 

12.8 

10.4 

4.0 

7.0 

Moderate-  Level 

10.2 

8.2 

5.0 

8.0 

Low-  Level 

8.0 

7.6 

6.0 

8.0 

Malta 

High-Levelt 
Moderate-  Level 
Low-  Level  t 

11.3 

10.8 

4.0 

5.0 

Norway 

9.3 

8.5 

3.0 

4.0 

*  Duration  is  the  actual  width  of  the  pulse  obtained  by  subtracting  the  interval  S  from 
the  definition  of  duration  in  Ref.  1.  c 

t  No  data  available. 


Table  I  illustrates  some  typical  results  using  the  above  criteria  which  will  be  of  use  to  us  later. 

It  is  of  some  interest  that  these  results  compare  favorably  with  visual  observations. 

It  should  be  noted  at  this  time  that  it  is  nowhere  suggested  that  the  simple  first-order  am¬ 
plitude  and  temporal  statistics  that  we  have  been  considering  completely  describe  ELF  noise. 
Nevertheless,  they  are  useful  in  that  any  proposed  model  must  at  least  result  in  statistics  which 
compare  favorably  with  these  easily  measured  quantities. 

m.  PROPOSED  MODEL 

An  intuitively  satisfying  model  for  impulsive  noise  processes  typical  of  the  ELF  environment 
consists  of  the  sum  of  a  low- density*  shot  process2'3  and  Gaussian  noise.  The  shot  process  ac¬ 
counts  for  the  randomly  occurring  high-level  pulses  while  the  additive  Gaussian  noise  represents 
the  low- level  background  noise.  A  convenient  representation  for  the  shot  process  is  that  of  a 
random  impulse  train  exciting  a  finite-state  linear  (usually  time  invariant)  system  as  first  pro¬ 
posed  by  Synder 4  Thus,  we  will  assume  that  the  noise  processes  of  interest  can  be  described 
by  the  state  equation 

4(t)  *  £*(t)  +bu(t)  ;  ,£(t0)  =jf0  (la) 

a*nd  the  output  relation 

y(t)  =  .  (lb) 

where  £  is  a  constant  n  x  n  matrix,  b  and  c  are  n- vectors  (specifically  column  vectors)  and 
<•  ,  •>  denotes  an  inner  product.  The  state  vector  x(t)  is  itself  an  n- vector  while  the  input  is  a 
scalar  of  the  form 

N(t) 

u(t)  =  £  ut  W* -  V  •  <2> 

i-i 

*  By  "low -density"  we  mean  that  the  average  Inter-arrival  time  between  exciting  Impulses  Is  greeter  than  any 
system  time  constants. 
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i.e.,  a  point  process.  Here  {N(t),  t  ^  0}  is  a  counting  process5  with  arrival  times  tj,  t2, . . . ,  and 
{uj}  is  a  sequence  of  independent  and  identically  distributed  (i.i.d.)  random  variables  with  common 
univariate  probability  density  function  (p.d.f.)  f(* ).  Finally,  the  observed  noise  is  given  by 

n(t)  =  y(t)  +  w(t)  ,  (3) 

where  w(t)  is  a  white*  zero-mean  Gaussian  noise  process  with  double-sided  noise  spectral 
density  =  Nq/2  watts/Hz  and  assumed  independent  of  the  driving  point  process  u(t). 

While  there  are  many  possible  state-space  realizations  of  the  system  described  by  Eq.  (1) 
leading  to  the  same  input-output  behavior,8  the  particular  realization  that  we  have  adopted 
(mainly  to  agree  with  Ref.  7)  is  such  that 


A  = 


-a. 


-a. 


-a 


n-1 


-  a_ 


(4) 


T  T 

while  b  =  (bt,  b2> , . . ,  bn>  and  finally  g  =  (c,  0,  0, . . . ,  0).  The  resulting  system  transfer  func¬ 
tion  of  the  linear  dynamical  system  generating  the  shot  process  is  easily  seen  to  be  given  by  a 
proper  rational  function  simply  related  to  the  coefficients  of  £  and  fe,  i.e., 

S)  _  y_  r_  t  a  i*ii 


(5) 


«s)  =<c,  [sI-A]'*b  > 


b.s"-1  +b2sn'2  +...b 

i  4  n 

sn +a1sn_1+...an 


where  I  is  the  n  x  n  identify  matrix  and  U(s)  and  Y(s)  are  Laplace  Transforms  in  the  complex 
variable  s  of  the  input  u(t)  and  output  y(t)  of  the  linear  dynamical  system.  The  resulting  state- 
space  realization  is  illustrated  in  Fig.  10.  An  alternative  realization  based  upon  the  system 
transfer  function  (Eq.  (5)]  which  will  prove  useful  in  later  digital  implementations  can  be  obtained 
by  observing8  that  the  matrix  fs|  -  in  Eq.  (5)  has  the  spectral  representation* 

°  SjL 


f-i-Af1  »  2  (s^Hs! 


k-1 


j 


(6) 


where  the  s^,  k  =  1,  2, . . .  n  are  eigenvalues  of  £  (poles  of  the  system  transfer  function)  and 


*  The  analysis  Is  easily  extended  to  colored  | 
t  We  have  assumed  that  the  poles  of  the  system  transfer  function  ate  distinct. 


mmv 


«,(*)■»!») 


Fig.  10.  Two  realizations  of  linear  dynamic  system,  (a)  State -space  realization; 

(b)  input/output  realization. 

2k  =  2H  §  Is!  -  Ar*ds  ;  k  =  1,  2, . . .  n  (7) 

J°k 

with  any  closed  contour  in  the  complex  plane  enclosing  only  the  pole  at  s^.  Substituting  this 
result  into  Eq.  (5),  we  obtain 


H(s)  =  £ 


"  A,. 


where 


mL  <£•  i-i-irV* 

=  ST  §r  ma)  ds 


is  simply  the  residue  of  H(s)  at  the  pole  sk>  Thus.  Eq.  (8)  is  merely  a  partial  fraction  expansion 
of  the  system  transfer  function. 

We  will  be  interested,  among  other  things,  in  the  amplitude  probability  distribution  (APD) 
of  the  process  n(t)  described  by  Eq.  (3).  To  this  end  consider  the  characteristic  function  (ch.f.) 
of  n(t)  given  by 

*n(v;  t)  A  E  {eivn(t>} 

.Efetv^.rft)^  B{#ivw(t)}  #  (10) 


where  we  have  made  use  of  the  assumed  independence  between  the  u(t)  and  w(t)  processes. 
Letting 

<yv)£E{elvw(t)} 


and 


Jv;t)—  E  {eiv<-*— ^>}  , 


(11) 


(12) 


we  have 


^n(v;t)  =  tpg(v)  ty(v;  t) 


(13) 


Note  that  ^(v)  is  merely  the  univariate  ch.f.  of  a  stationary  zero-mean  Gaussian  process  so  that 


Vg(v)  =  e 


-K2v2 


(14) 


where  has  been  described  previously  and  is  independent  of  t.  To  evaluate  ip y(v,  t),  note  that 
the  solution  to  (la)  is  given  by^ 


x(t)  =  *(t.  tQ)  x(tQ)  + 


5"  - 


bu(r)  dT  , 


(15) 


where  * (t,  t)  =  exp  {  A(t  -  r)}  is  the  state  transition  function.  From  Eq.  (2)  this  last  expression 
can  be  written  as 


N(t) 


x(t)  =  *<t,  tQ)  x(tQ)  +  2)  V  feuk 


(16) 


k=N(tQ)+l 


so  that  finally 


’PJv.t)  =  E  {exp[iv<c,  |(t,t  )  x(t  )>J}  E  jexpfiv  2  <£,*<t,  V  -^"k  II  *  <17> 

'  «•  k=N(t0)+l  ■»' 

Let  us  give  consideration  to  the  evaluation  of  the  second  of  these  expectations.  This  quantity  is 

easily  seen  to  be  the  ch.f.  of  a  random  sum  of  independent  although  not  necessarily  Identically 

distributed  random  variables.  If  the  counting  process  {N(t),  t  >0}  is  assumed  Poisson4'  with 

.5 


intensity  X,  it  can  easily  be  shown  that 
N(t) 


jexp  |iv  2  <£.  *(*•  V  ukj  |  =  e*P  r)  -  ij  dr  j  , 

k  N(tQ)+l  o 


(18) 


where 


S»(v;  r)  *  E  {expflv^c ,  *(t,  r)  b>  u)  } 


(19) 


*  It  U  of  tome  Interest  to  extend  these  results  to  the  ease  of  an  aibltrary  renewal  counting  process. 
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is  the  ch.f.  of  the  random  term  appearing  in  the  summation  of  (14).  Here  the  expectation  is 
with  respect  to  the  only  random  quantity  u  so  that 

<p(v;  t)  =  ^u(v<c ,  £(t,  r)  b>)  (20) 

with  Vu(  •  )  the  univariate  ch.f.  of  the  random  amplitudes  {u^}  possessing  common  p.d.f.  f(*  ). 

It  will  be  convenient  to  let 

h(t,  r)  =  <c,  *(t,  t)  b>  ,  (jD 

noting  that  h(t,  r)  =  h(t  -  r)  due  to  the  time  invariance  of  *(t,  r),  Indeed,  h(t)  is  easily  seen  to  be 
the  inverse  Laplace  Transform  of  the  system  transfer  function  H(s)  in  Eq.  (5)  and  thus  represents 
the  impulse  response  of  the  linear  time-invariant  system. 

Similarly,  for  the  first  expectation  in  (17)  we  obtain 

E  {exp(iv<c,£(t,to)  x(tQ)>}  =  E  (exp[i<v*t{t.t0)  c,x(tQ)>]  }  £  (WNt.t  )  c)  (22) 

—  o' 

with  ^x(t  )( *  *  the  multidimensional  ch.f.  of  the  initial  state  vector  x(tQ),  and  where  #*(t,t  )  de¬ 
notes  the0 conjugate  transpose  of  the  state-transition  matrix.  Thus,  using  the  preceding  results 
in  (17),  we  obtain 

^y(v;  t)  =  *x(y  (v*f(t,  tQ)c)  exp  |x  £  [*u<vh(t,  r))  -  1)  drj  (23) 

for  the  ch.f.  of  the  shot  component.  For  a  specified  distribution  of  the  initial  state  vector  x(t  ) 
and  the  p.d.f.  f(  • )  of  the  random  amplitude,  the  ch.f.'s  *  ,.  .(  • )  and  tpjl  ■ )  are  determined.  The 
evaluation  of  tfy(v;  t)  is  then  relatively  straightforward,  although  possibly  tedious,  for  specified 
system  dynamics.  The  quantity  tfy(v;  t)  has  been  evaluated  for  some  special  cases  in  a  previous 
memorandum.9 

Also  of  some  interest  is  the  joint  ch.f.  of  the  shot  component  at  times  tj,  t2  with  t2  >  tj 
say.  By  a  derivation  similar  to  the  above  it  is  easily  shown  that  this  quantity  is  given  by 

V  V  "  *x(tQ)  <^1** <V  V  +  v2*t(t2.  tQ)j  c)  exp  |x  £*  I*u{vlh(t,  r) 

"  o 

+  v2h(t2.  r)}  -  1]  drj  exp  |x  [*u(v2h(t2.  r»  -  lj  drj  ,  (24) 

.where  y  =  (Vj,  v2).  This  expression  can  again  be  evaluated  for  specified  distributions  on  x(to) 
and  u  and  assigned  system  dynamics. 

It  will  be  convenient  in  what  follows  to  assume  that  we  will  be  interested  only  in  times  that 
are  sufficiently  far  removed  from  tQ  so  that  initial  transient  disturbances  have  been  dissipated. 

For  all  such  times  t  we  have  approximately  **(t,  tQ)  c  "  0,  the  null  vector  so  that  the  following 
approximations  are  valid 

*y(v;t)  *  exp  jx  f*u(vh({))  -  1)  d(j  (25) 


and 


<yv;  tr  t?)  *  exp  |A  f^u(v1h({)  +  v2h({  +  r)}  -  i} dfj 

X  exp  |a  JJ  f^(v2h({))  -  1 J  d{  J  ,  (26) 

where  r  =  tj.  With  y (t)  -  ^c,*(t))  representing  the  shot  component,  it  follows  that 

E  (y(t)}  =  AE  {u}  f  h({)  df  .  (27) 

J0 

Similarly,  we  have 

_BO 

var{y(t)}  =  AE  {u2}  \  h2(f)  d{  (28) 

and  for  the  autocorrelation  function 

R( r)  =  E  {y(t)  y(t  +  r)}  •  A2E2  {u}  [  j*  h«)  d(J  2  +  AE  {u2>  J  h({)  h(«  +  r)  df  (29) 

with  corresponding  spectral  density  function 

S(w)  =  2xA2E2  {u}  H*«»  Hu)  +  AE  {u2}  |H(jw)!2  ,  (30) 

where  H(ju)  is  the  Fourier  Transform  of  h(t),  i.e.,  the  system  transfer  function  H(s)  evaluated 
at  s  =  ju  provided  that  the  imaginary  axis  is  within  the  region  of  absolute  convergence  of  Us). 
Table  n  summarizes  the  properties  of  several  linear  systems  which  will  be  of  interest  to  ua  in 


TABLE  II 

MEAN  AND  VARIANCE  OF  SHOT  COMPONENT 
FOR  SOME  USEFUL  SYSTEM  TRANSFER  FUNCTIONS 


H(e) 

h(t) 

E  (y<t)} 

Var{y<t)} 

c 

m 

ce'^u.^t) 

AcE  {u} 

a 

Ac2E  {u2> 
- To 

c 

—  s{h  /,if  «| 

AcE  {u> 

Ac2E  fu2} 

(a  +  a)*  +  ul 

* sin  w X  U  4**1 

Cl)  “» 

az  +  «z 

4a<az  +  wZ) 

c(s  +  a) 

ce~at  cosut  u_^(t) 

AacE  (u) 

Ac2«2E  (u2) 

(a  +a)*  +w2 

az  +«z 

4a(a*  +  «z) 

what  follows.  We  will  for  the  most  part  be  interested  only  in  the  ease  where  S  {u)  *  0.  From 
Table  II,  then,  it  is  a  simple  matter  to  choose  the  value  c  for  specified  system  dynamics  to 
achieve  a  given  value  of  y  defined  by 

yl  .  (34) 

i.e.,  y  represents  the  ratio  of  the  rms  value  of  the  shot  noise  component  to  that  of  the  additive 
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background  Gaussian  noise  and  will  prove  to  be  an  important  parameter  of  the  proposed  model. 
In  summary,  then,  the  parameters  of  the  proposed  model  which  must  be  specified  are: 


(a)  The  p.d.f.  f(  • )  of  the  pulse  amplitudes, 

(b)  The  intensity  X  of  the  Poisson  point  process, 

(c)  The  linear  system  dynamics, 

(d)  The  ratio  y  of  the  rms  values  of  the  shot  and  Gaussian  components. 

2 

Observe  that  the  choice  of  background  noise  level  <r  is  somewhat  arbitrary  and  merely  amounts 
to  a  trivial  rescaling  of  the  process  n(t).  As  a  result,  we  will  assume  a*  =  i  in  what  follows. 

In  the  following  sections  we  will  discuss  the  choice  of  the  above  parameters  to  match  observed 
ELF  noise. 

IV.  DIGITAL  REALIZATION  OF  MODEL  FOR  ELF  H- FIELD  NOISE 

We  will  be  particularly  concerned  with  fitting  the  model  described  in  the  preceding  section 
to  observed  H-field  data  and  obtaining  a  suitable  digital  realization  of  this  model.  The  situation 
is  complicated  somewhat  by  the  fact  that  the  available  noise  data  have  been  recorded  with  a  loop 
antenna,  thus  measuring  the  time  derivative  of  the  ambient  H-field.  To  compensate  for  this, 
the  raw  data  are  processed  by  a  filter  with  an  attenuation  characteristic  exhibiting  a  6  dB/octave 
roll-off  out  to  350  Hz.  This  approximate  integration  reconstructs  a  facsimile  of  the  actual 
H-field  data.  We  will  approximate  this  filter  by  a  moving  average  operation  performed  on  the 


u(t) 

LINEAR  DYNAMICAL 

SYSTEM 

i  /'  H. 

*  va 


Fffl.  11 •  Modal  of  obforvod  H-flold  do  to. 

n(t)  process  as  illustrated  in  Fig.  11,  The  resulting  process  nit)  is  then  related  to  the  n(t) 
process  of  the  previous  section  by 


n(f)  dr  . 


This  operation  can  be  considered  as  the  convolution  of  the  input  n(t)  with  the  Impulse  response 

hjlt)  •  £  fu.jlt)  -  u.jlt  -  «)J  .  (33) 

wrtere  u_j(t)  is  the  unit  step  function.  The  corresponding  squared- magnitude  response  Is 
given  by 

'  <«> 

which  is  down  3  dB  at  u  «  «  &  2.78/s.  ft  is  reasonabls  then  to  adjust  the  average  interval  <  to 
result  in  a  "bandwidth"  «Q  *  2t  x  350  rad /sec.*  We  will  assume  then  that 

*  Later  studies  have  Indicated  that  ft  should  be  made  somoehet  larger  than  this  value  to  result  In  spectral 
characteristics  which  mare  closely  resemble  observed  ELF  H -field  noise. 


WwiWsse,.* 


(35) 


6  = 


JL 

2* 


<w> 


in  what  follows. 

The  digital  simulation  of  the  input/output  behavior  of  the  linear  dynamical  system  is  deter¬ 
mined  from  the  following  considerations.  Observe  that  the  solution  to  the  state  equation  (la)  is 
given  for  any  t  >  r  by 


x(t)  =  *(t,  t)  x(  r)  +  J  £(t,  {)  bu({  )  d{ 


(36) 


where  *(t,  r)  =  exp{A(t  -  r)}  as  has  been  described  previously.  Thus,  if  we  are  interested 
only  in  the  discrete  times  t  =  kA,  k  =  0, 1,  2, ... ,  where  A  is  the  sampling  interval,  we  have 
upon  letting  x^  =  x(kA)  the  discrete  recurrence  relation 

-k+i  =--k +-kuk+l  :  -o  =  -  ’  (37a> 


yk  *  <£»*k>  .  (37b) 

where  F  =  exp  (A  A},  and 

0  if  no  impulse  occurred  in  the  interval  (kA,  (k  +  i) 

Uj  if  the  i^  impulse  occurred  at  time  t^,  kA  <  tj  ^  (k 

Finally,  b  is  defined  only  for  nonzero  u^  values  as  b  k  =  exp  {A[(k  +  1) 
z-transforms  in  (37),  we  have  for  the  discrete  transfer  function 

H(z)  =  =  <c.  [I- z^F]"1  b>  .  (39) 

U(z)  “ 

where  U(z)  and  Y(z)  are  the  z~transforms  of  the  input  and  output  sequences  {u^}  and  {yk}. 
respectively.  We  shall  assume  that  the  sampling  Interval  A  is  small  relative  to  any  time  con¬ 
stants  associated  with  the  linear  system  so  that*  b  »  b.  Again  making  use  of  the  spectral 
representation8  we  have 


A) 

(38) 

+  1)  A  . 

A  -  tj]}  b.  Thus,  taking 


z 


k=l 


(40) 


where  the  V  k  =  4»  2'  •  •  *  n  are  again  poles  of  H(s),  and  the  £](  are  defined  as  in  Eq.  (6)  and  the 
sequel.  It  follows  that  the  discrete  transfer  function  H(s)  can  be  expressed  as 


»  A. 

H(z)  -  Z  - S-j-X  (41) 

-1 

k«l  1  -  z  1  e  K 

with  Ak  k  *  1,  2, . . .  n,  again  the  residues  of  H(s)  at  the  distinct  poles  s^,  k»i,2,.,.u,  This 
is  the  essence  of  the  impulse  Invariance  technique  of  digital  filtering10  and  results  in  a  parallel 
realization  of  the  discrete  system  described  by  ( 37).  In  the  present  study,  we  have  concentrated 
upon  the  three  linear  systems  described  in  Table  II  for  which  the  parallel  digital  realizations 


•  At  any  rate.  In  the  digital  simulation  the  Impulses  occurring  at  kA  <  tj  £  (k  +  I)  A  will  always  be  placed  at  the 
right-hand  boundary  of  such  Intervals  so  that  this  apprcKlmatlon  Is  exact. 
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corresponding  to  (41)  are  readily  obtained.  A  FORTRAN  program  for  a  complete  digital  simula¬ 
tion  of  the  proposed  model  has  been  written  and  used  extensively  in  this  study. 

V.  CHOICE  OF  PARAMETER  VALUES 


The  first  quantity  to  be  specified  for  the  proposed  model  is,  of  course,  the  p.d.f.  f(  • ) 
of  the  individual  pulse  amplitudes.  In  principle,  this  information  could  be  obtained  from  the 
results  of  Section  III  by  computing  the  APD  from  the  derived  expressions  for  the  ch.f.  for  a 
specific  choice  of  f(  • )  and  comparing  the  results  with  measured  APD's.  The  difficulty  in  actually 
computing  the  APD  except  in  special  cases,  however,  precludes  this  approach.  As  an  alterna¬ 
tive,  let  us  note  that  the  APD  at  high  levels  should  essentially  be  determined  by  the  large  noise 
peaks  and  hence  should  resemble  f(*  ).  Let  us  then  choose  a  p.d.f.  whose  amplitude  distribution 
at  large  levels  is  similar  to  that  observed  in  measured  APD's.  We  have  considered  two  such 
zero-mean  distributions  which  appear  to  satisfy  this  requirement:  the  double-sided  power- 
Rayleigh  distribution  defined  by 


f(u)  = 


0  <  a  <  2 


(42) 


and  the  double-sided  log-normal  distribution  defined  by* 

»2 


f(u) 


2  ^/2^^ a  lu  i 


exp 


_  (1 n  Ju  |  -  nr 

Us 


(43) 


2 

where  p  and  a  are  the  mean  and  variance,  respectively,  of  a  Gaussian  variate  g  for  which 
u  =  e®.  In  particular,  we  will  assume  that  p  =  0  in  what  follows.  The  power- Rayleigh  distribution 
has  a  variance  R  *1X1+  2/a),  where  It  • )  is  the  Gamma  function  while  the  log-normal  distribu¬ 
tion  possesses  a  variance  exp  {2a2}.  Thus,  the  power- Rayleigh  distribution  is  completely  defined 
by  the  scale  parameter  RQ  and  the  exponent  «.  A  plot  of  the  resulting  APD  for  the  power- 
Rayleigh  distribution  is  illustrated  in  Fig.  12  for  several  values  of  a.  By  comparison  with  the 
high-level  region  of  the  measured  APD  plots  given  in  Figs.  4  through  6,  values  cf  cy  in  the  range 
0.25-$  a  4  0-  50  have  been  found  appropriate.  Similarly,  we  have  found  it  useful  to  characterize 
the  log-normal  distribution  by  the  scale  parameter  a  and  the  related  quantity 


Vd  = 10  lo*i0 


.JLfo2? , 

E2  {  |  u  | } 


(44) 


=  4. 34  a2  . 


Plots  of  the  resulting  APD  for  several  values  of  Vd  are  illustrated  in  Fig,  1 3.  Here  it  appears 
that  values  of  in  the  range  10  Vd  15  prove  sufficient. 

Now  let  us  consider  appropriate  choices  for  the  Intensity  X  of  the  point  process.  This  in¬ 
formation  can  be  obtained  directly  from  the  data  of  Table  I.  For  definiteness  we  will  take  X  as 
the  average  number  of  high-level  pulses  per  second  from  Table  I  with  Sc  «5  msec  as  this  choice 
appears  to  correlate  most  favorably  with  visual  observations.  Additional  information  on  system 


•Each  of  those  distributions  eon  be  obtained  by  forming  the  random  variable  u  «*  pr,  where  r  has  a  conventional 
one-sided  distribution  of  either  type  and  p  is  a  random  variable  taking  values  el  with  equal  probability. 
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Fig.  12.  APO  of  double-sided  power- Rayleigh  distribution. 


Fig.  13.  APO  of  doubts -sided  log-normal  distribution.  1 
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dynamics  can  be  gleaned  from  Table  I.  In  particular,  each  of  the  linear  system  transfer  func¬ 
tions  in  Table  II  possesses  an  impulse  response  with  a  well-defined  duration  (or  time-constant) 

Tc  =  i/a  seconds.  It  then  proves  convenient  in  each  case  to  adjust  the  parameter  a  to  provide 
an  average  pulse  duration  as  given  by  the  measured  data  in  Tab)e  L  Again  for  definiteness  we 
take  the  average  pulse  duration  with  Sc  =  5  msec  as  this  choice  correlates  most  favorably  with 
visual  observations.  For  system  transfer  functions  with  complex  poles,  it  has  been  determined 
that  w~  2000  rad/sec  should  prove  sufficient  in  each  case  considered.  Thus,  all  that  remains  is 
to  develop  an  appropriate  rationale  for  choice  of  the  quantity  y  given  by  (31).  This  in  turn  fol¬ 
lows  easily  from  consideration  of  the  effect  of  the  parameter  y  on  a  typical  APD  plot.  As  we 
have  previously  observed,  the  APD  is  controlled  mainly  by  the  background  Gaussian  noise  in  the 
small  signal  regime.  Thus,  by  observing  the  amount  by  which  the  straight-line  segments  of 
the  APD  plots  in  the  small-signal  regime  are  depressed  below  the  APD  corresponding  to  Gaussian 
noise  alone  (Fig.  7),  we  should  be  able  to  obtain  a  rough  estimate  of  y.  In  particular,  observe 


2  2  2 

where  we  have  made  use  of  the  fact  a  *  a  (1  +  >•  ),  Thus,  the  amount  by  which  the  APD  should 
be  depressed  below  that  for  Gaussian  noise  alone  is  given  by 


AdB  =  10  log1Q(l  +y‘)  .  (46) 

It  should  be  noted  that  this  approximation  can  be  expected  to  hold  only  for  y  sufficiently  greater 
than  unity.  As  an  example,  comparing  Fig.  4  with  Fig.  7,  it  is  observed  that  the  low-level 
Gaussian  segment  of  the  former  is  depressed  by  AdB  «  13  so  that  y  =  4.4,  In  cases  where  the 
approximation  fails,  resort  has  been  made  to  visual  observation  of  the  resulting  APD  plots  for 
choice  of  y. 

Table  III  indicates  the  appropriate  choice  of  parameter  values  for  each  of  the  models  con¬ 
sidered  and  has  been  generated  by  the  above  procedure.  As  an  example  of  the  type  of  results  to 
be  obtained,  consider  the  high-level  Saipan  data.  Figures  14  through  19  illustrate  typical  noise 
waveforms  obtained  from  the  proposed  model  with  parameter  values  taken  from  Table  IIL  These 
waveforms  should  be  compared  with  that  in  Fig.  1.  Observe  that  the  most  favorable  visual  corre¬ 
spondence  is  obtained  with  a  linear  system  possessing  a  pair  of  complex  conjugate  poles.  We 
have,  nevertheless,  decided  not  to  neglect  entirely  the  single-pole  model  as  it  is  felt  that  this 
simplified  model  may  have  some  utility  in  later  studies  -  in  particular,  in  the  digital  computer 
simulation  of  optimum  detector  structures.  It  has  been  determined  that  the  resulting  APD  plots 
for  the  artificially  generated  ELF  noise  do  not  exhibit  noticeable  differences  for  the  different 
linear  system  used  to  generate  the  shot  process,  provided  that  the  other  parameters  remain 
fixed.  Figures  20(a)  and  20(b)  illustrate  the  correspondence  between  the  APD's  of  artificially 
generated  and  recorded  high-level  Saipnn  data  when  the  parameters  of  the  proposed  model  are 
cho«.on  from  Table  III.  As  can  be  observed,  the  correspondence  is  extremely  good.  Similar 
correspondences  have  been  obtained  with  the  Malta  and  Norway  data. 

As  a  final  comment,  it  should  be  noted  that  some  estimate  of  the  processing  gains  to  be 
realized  with  typical  nonlinear  receiver  structures  operating  upon  noise  as  described  by  this 
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Fig.  14.  Simulated  noise  waveform:  single-pole  model  with  power-Rayleigh 
distributed  pulse  amplitudes. 


Fig.  15.  Simulated  noise  waveform:  two -pole  modal  with  power-Rayleigh 
distributed  pulse  amplitudes. 
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Fig.  16.  Simulated  noise  waveform:  two-pole,  single  xero  model  with  power-Rayleigh 
distributed  pulse  amplitudes. 


Fig.  17.  Simulated  noise  waveform i  single-pole  model  with  log-normal 
distributed  pulse  amplitudes. 
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Fig.  20.  Corratpondanea  between  APO  of  ELF  nolle  modal  and  maaiured  data 


model  can  be  obtained  from  the  data  of  Table  III.  In  particular,  consider  a  detector  of  the  hole- 
punching  variety  with  threshold  set  just  above  the  background  Gaussian  noise  level*  The  fraction 
of  signal  energy  lost  in  any  1 -second  interval  is  approximately  XTp  which,  from  the  data  of 
Table  III,  is  negligible  in  each  case.  The  reduction  of  noise  power,  on  the  other  hand,  is  approx¬ 
imately  given  by  A^g  as  defined  in  (46).  Thus,  the  processing  gain  is  approximately  A^g  which, 
in  each  of  the  cases  considered,  has  been  included  in  Table  III.  This  estimate  comes  surpris¬ 
ingly  close  to  the  processing  gain  measured  with  corresponding  recorded  noise  data. 

VI.  SUMMARY  AND  CONCLUSIONS 

A  model  for  ELF  noise  has  been  proposed  and  studied  in  detail.  A  procedure  has  been  given 
for  fitting  the  free  parameters  of  this  model  to  observed  noise  data.  It  is  felt  that  the  model 
satisfies  all  the  diserata  listed  in  Section  L  Further  studies  will  be  concerned  with  the  determi¬ 
nation  of  optimum/suboptimum  detector  structures  for  reception  in  noise  as  described  by  this 
model  and  the  evaluation  of  the  resulting  performance. 

There  are  several  directions  in  which  the  proposed  model  could  be  generalized  to  more 
closely  represent  observed  ELF  noise.  For  instance,  the  Poisson  point  process  could  be  re¬ 
placed  by  one  which  would  result  in  some  correlation  between  noise  burst  intervals  as  has  been 
obse,  zed  in  recorded  data.  Another  possibility  is  to  consider  linear  systems  with  multiple- 
order  poles  for  generating  the  low-density  shot  processes.  This  would  result  in  pulse  waveforms 
possessing  observable  precursors  which  again  have  been  observed  in  recorded  data.  Such  re¬ 
finements  were  not  undertaken  in  the  present  study. 
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